The invariant joint distribution of a stationary random field and its derivatives: 
Euler characteristic and critical point counts in 2 and 3D 
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The full moments expansion of the joint probability distribution of an isotropic random field, its 
gradient and invariants of the Hessian is presented in 2 and 3D. It allows for explicit expression for 
the Euler characteristic in ND and computation of extrema counts as functions of the excursion set 
threshold and the spectral parameter, as illustrated on model examples. 
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Random fields are ubiquitous phenomena in physics 
appearing in areas from turbulence to the landscape of 
string theories. In cosmology, the large scale distribution 
of matter (LSS) and the sky-maps of the polarized Cos- 
mic Microwave Background (CMB) radiation - two focal 
topics of current research - are described as, respectively, 
3D and 2D random fields. Modern view of the Universe, 
developed primarily through statistical analysis of these 
fields, points to a Universe that is statistically homoge- 
neous and isotropic with a hierarchy of structures arising 
from small Gaussian fluctuations of quantum origin. 

While the Gaussian limit provides the fundamental 
starting point in the study of random fields [1-3J, non- 
Gaussian features of the CMB and LSS fields are of great 
interest. CMB inherits high level of Gaussianity from the 
initial fluctuations, and small non-Gaussian deviations 
may provide a unique window into the details of processes 
in the early Universe. The gravitational instability that 
nonlinearly maps the initial Gaussian inhomogeneities in 
matter density into the LSS, on the other hand, induces 
strong non-Gaussian features culminating in the forma- 
tion of collapsed, self-gravitating objects such as galaxies 
and clusters of galaxies. At supercluster scales where 
non-linearity is mild, the non-Gaussianity of the mat- 
ter density is also mild, but still essential for quantita- 
tive understanding of the filamentary Cosmic Web [3] 
in-between the galaxy clusters. 

The search for the best methods analyzing non- 
Gaussian random fields both in weak and strong regimes 
is ongoing. We focus on the statistics of geometrical and 
topological properties of the field that includes the Eu- 
ler characteristic of excursion sets [SHE] (and the rest of 
Minkowski functionals [5]), the density of extremal points 
and statistics of critical lines, or the skeleton [10l[TT]. In 
this paper we present the formalism for computing such 
geometrical statistics for the homogeneous and isotropic 
mildly non-Gaussian fields that can be represented or 
approximated by a Gram-Charlier expansion around the 
Gaussian limit, to an arbitrary order of the expansion. 

A statistically homogeneous ND random field x is 
fully characterized by the joint one-point distribution 
function (JPDF) of its value and its derivatives P(x), 
x = (x, Xi, Xij, Xijk, ■ ■ •)■ We consider non-Gaussian fields 
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around the Gaussian G(x) = (2tt)- jv / 2 |C|- 1 / 2 exp(-±x- 
C _1 • x) that is arranged to match the mean and the 
covariance C = (x ® x) of the x variables. We con- 
sider all the variables to be defined as having zero mean 
and unit variance, Cu — 1. The correction to the 
Gaussian approximation is the series in Hermite tensors 
h„(x) = (-l) n G- 1 (n)d n G(x;)/dx n of rank n with coeffi- 
cients constructed from the moments, (x n ) GC = (h„(x)). 

The statistics of Euler characteristic and extremal 
points [T] , and the basic description of the critical lines of 
the field in the stiff approximation [llj , require to know 
the JPDF only up to the second derivatives. In the Gaus- 
sian limit, in this case, the only non-trivial covariance 
parameter is the cross-correlation between the field and 
the trace of the Hessian 7 = — (xTi(xij)} [HUH]- 

The JPDF in the form of Eq. is not ideal to study 
critical sets statistics since the coordinate representation 
masks the isotropic nature of the statistical descriptors. 
The pioneering works [51 [7] , where the first correction to 
Euler characteristic was computed, demonstrate the aris- 
ing complexities. Instead we develop the equivalent of the 
Gram-Charlier expansion for the JPDF of the field vari- 
ables that are invariant under coordinate rotation. Such 
distribution can be computed via explicit integration of 
the series Eq. Q over rotations, however we obtain it di- 
rectly from general principles: the moment expansion of 
the non-Gaussian JPDF corresponds to the expansion in 
the set of polynomials which are orthogonal with respect 
to the weight provided by the JPDF in the Gaussian limit. 
Thus, the problem is reduced to finding such polynomials 
for a suitable set of invariant variables. 

The rotational invariants that are present in the prob- 
lem are: the field value x itself, the modulus of its gra- 
dient, q 2 = Xi 2 and the invariants of the matrix of 
the second derivatives x^-. A rank N symmetric matrix 
has N invariants with respect to rotations. The eigen- 
values Xi provide one such representation of invariants, 
however they are complex algebraic functions of the ma- 
trix components. An alternative representation is a set of 
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invariants that are polynomial in Xij, with one indepen- 
dent invariant polynomial per order, from one to N. A 
familiar example is the set of coefficients, I s , of the char- 
acteristic equation for the eigenvalues, where the linear 
invariant is the trace, I± — ^ZjAj, the quadratic one is 
I2 = Ylj^j AjAj and the N-th order invariant is the deter- 
minant of the matrix. In — J\ i A;. Aiming at simplifying 
the JPDF in the Gaussian limit |T7| - (e.g., [Sin], Ap- 
pendix A) we use their linear combinations J s 
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where J s >2 are (renormalized) coefficients of the charac- 



terstic equation of the traceless part of the Hessian and 
are independent in the Gaussian limit on the trace J\. 

Let us consider the 2D and 3D cases explicitly. Intro- 
ducing Q — (x + 7 Ji)/ ^\ — 7 2 in place of the field value 
x we find that the 2D Gaussian JPDF G 2D (C, <? 2 , Jl, J2), 
normalized over d(dq 2 dJidJ2, has a fully factorized form 



G 2D = x- exp 
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Used as a kernel for the polynomial expansion, G 2 d leads 
to a non-Gaussian rotation invariant JPDF in the form 
of the direct series in the products of Hermite, for £ and 
Ji, and Laguerre, for q 2 and J2, polynomials: 



P2D((,q 2 ,Jl,J2) = G 2D 



00 i+2j+k+2l=n 

E E 

n— 3 iJ,k,l—0 



i\ jl kl 11 



/-i„23 t k t I 

C, q J\ J2 



H{() Ljfq^HkiJjL^Jz) 



(4) 



where ^Y^tj'k ~i=$ 2l ~ n stands for summation over all com- 
binations of non- negative i,j,k,l such that i + 2j + k + 21 
adds to the order of the expansion term n. The coeffi- 
cients of the expansion are combinations of moments 
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that vanish in the Gaussian limit. In the lowest n = 3 
order they coincide with the cumulants of our variables. 

In ND the Gaussian JPDF, Gnd(0 q 2 , Jl, Js>2), re- 
tains complete factorization with respect to C, q 2 , Ji 
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so in the moment expansion this sector always gives 
rise to Hermite Hi((), H k (Ji) and generalized Laguerre 
j^(n 2)/2^^2 polynomials. However, the distribu- 
tions of the rest J s > 2 contained in Q{J S >2) are coupled. 



Specifically, in 3D, Q( J 2 , J 3 
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is distributed uniformly between — J 2 and J 2 . Let us 
denote the orthogonal polynomials in these two variables, 
Fim{ J2, J3), where I is the power of J 2 and m is the power 
of J3. They obey the orthonormality condition 
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We shall not give the full theory of these polynomials 
here, but note that one can construct them by a Gram- 
Schmidt orthogonalization procedure to any given order. 
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J3 are sufficient to obtain the general expression for 
the Euler characteristic of the excursion sets of the field 
to arbitrary order, and to calculate the critical point and 
skeleton statistics to quartic order. 

Hence, we write the moment expansion for invariant 
non-Gaussian JPDF, P 3 u{(,q 2 , Ji, J2, J3) as a series in 
the power order of the field, n in the form 
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The first high order term contains J2 but not J3, the 
second is linear in J3 and contains no J2, while for where 
we left the normalization coefficient q m undetermined, 
contains all the remaining combinations of J 2 and J3. 
In this last term, C01 = 0, so that the first contribution 
to it, I = m = 1, is of the fifth power of the field oc 
J2J3. The moment combinations that give the expansion 
coefficients are found analogously to 2D case. 

The formulas Q and (|6| provide the joint probability 
function of a field and its derivatives up to second order 
in terms of the invariant variables to an arbitrary order 
in the moment expansion. They allow to easily compute 
any rotation-invariant statistics that depend exclusively 
on these descriptors of the field. In particular, one can 
compute the Euler characteristic of the excursion sets of 
the field analytically, the density of the extrema, and the 
properties of the critical lines that describe the skeleton of 
the field, as we describe in detail in the follow-up paper. 

Determining the topological Euler characteristic, Xi 
and the density of the extrema in the high excursion 
sets x > v of the random field are two related classi- 
cal problems in the study of the geometry of a random 
field. Both are reduced [T] to the evaluation under the 



extremal condition q 2 = of the statistical average of 
the Gaussian curvature of the field surface given by the 
invariant the difference being the various conditions 
set on the signs of the eigenvalues of the Hessian. For 
the Euler characteristic integration is unconstrained 

^ = (-1) N Jjx jdq 2 q N -H N D {f) jf[dJ s P ND I N , 

while to find the number density of extremal points of 
different types one integrates over the regions in J s space 
where particular signature of the eigenvalues is main- 
tained. Since I s are just low order polynomials of the 
J s variables, Eqs Q and (|6| are well suited to perform 
the required calcualtions. 

The Euler characteristic can be computed completely 
by noting from Eq. ^ that 7jy depends only linearly on 
J s > 2 (e.g., h = 4f ( J i 3 - 3J i J 2 + 2 ^3) in 3D), hence all 
terms in JPDF of higher order in J s >2 do not contribute. 
The 2D and 3D results [18. can be combined in a very 
compact form if one re-expresses the coefficients back in 
terms of the field x itself and the invariants I s 
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where i — 0, s — N terms have been combined into the 
boundary term cx \{~ 00 ) fixed by the topology of the 
manifold and should be omitted from the sum. Iq = 1. 

Calculation of the extrema number density produces 
an analytical result only in 2D. Even then the result to 
fourth order is already too complicated to be reproduced 
here and is deferred to |13j . Instead here we demonstrate 
our results on two model non-Gaussian fields. 

For the a 2D case we generate a Gaussian random field 
xg with the scale-invariant power spectrum cx fc~ 1,5 , 
smoothed with a Gaussian filter of 5-pixel width. The 
underlying grid is then displaced by aWA~ 1 xa and the 
initial field is resampled on this grid, producing a one- 
parameter non-Gaussian field, x^- a K This toy model 
is able to produce high level of non-Gaussianity with- 
out generating high field excursions, as demonstrated in 
Fig. |l] by a 2048 2 realization of the field with a = 0.5. 

In 3D, we use scale-invariant (oc simulations of 

cosmological density evolved to the mildly non-linear 
stage of gravitational instability where the variance of 
inhomogeneities is a — 0.1 of the mean density. 

In Fig. [2] we plot the extrema count [IH] and the Euler 




FIG. 1: A section of strongly non-Gaussian 2D field. 



characteristic computed within our formalism to 3rd and 
4th order respectively, using the moments measured nu- 
merically in our toy simulations. The left column demon- 
strates the avenues our formalism opens for the theoret- 
ical study of critical points of non-Gaussian fields. The 
right column shows, on the example of the Euler charac- 
teristic, the importance of accounting for high moments 
of the expansion for some models of non-Gaussianity, as 
well as the convergence process of the expansion for the 
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FIG. 2: The number of extrema (left) and the Euler character- 
istic (right) for the 2D toy model (top) and 3D gravitational 
collapse model (bottom). In dimensional units both quanti- 
ties are given per volume defined as in [TT]- The shaded 
bands correspond to 2a variations over the mean measure- 
ments of 100 (2D) or 25 (3D) realizations, while the curves 
give the 3 r and 4 th order predictions. A higher order cor- 
rection is required to fit the strong non-Gaussian field shown 
in Fig. [I] In 2D only the correction to the Gaussian limit is 
shown. In 3D, the Gaussian prediction is shown as a thin line. 

gravitational instability models. 

The ability to analyze the full series of the moment ex- 
pansion for geometrical statistics is important from sev- 



eral points of view. In cosmological applications, e.g. a 
mildly non-linear gravitational instability as simulated 
in 3D in Figure [2j or lensing of CMB maps, one deals 
with perturbative deviation from Gaussianity governed 
by the small value of the variance a. In this case Gram- 
Charlier series can be rearranged in (asymptotic) Edge- 
worth 14 power series in a |20j and the full expansion 
can be used to obtain an increasingly accurate descrip- 
tion of the statistics, and conduct estimates of residual 
errors. Outside of established cosmological perturbation 
theory, theoretical formalism that deals globally with all 
the terms is all the more important. It may allow us 
to design a customized truncation criteria, which may 
not be uniform across all the statistics, or even for all 
the partial contributions to a given statistical descriptor 

In the cosmological context, non- Gaussianity of extra- 
galactic fields can be used to constrain the dark energy 
equation of state via 3D galactic surveys, or shed light 
on the physics of the early Universe through 2D CMB 
maps. The predictions of theoretical models are often 
given through the hierarchy of the differences between 
the moments to their Gaussian limit. Higher order mo- 
ments are generally difficult to test directly in real-life 
observations, due to their sensitivity to very rare events. 
The geometrical analysis of the critical events in the field 
provides more robust measures of non-Gaussianity and 
is becoming an active field of investigation [T51 [TO] . The 
present paper advances the formal development of this 
subject, opening ways to expand the list and sharpen 
statistical tools that may detect unique signatures of fun- 
damental processes in our Universe. 
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